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ABSTRACT 

This study consists of the development of a computer program to 
numerically solve the space and energy dependent multigroup neutron 
diffusion equations in a bare homogeneous fast reactor core or reactor 
material assembly. 

The resulting program is unique in that it was designed for future 

use by Naval Postgraduate School students undertaking experimental 

studies in neutron diffusion with limited time to determine numerical 

solutions for verification of their results. 

The equations are solved iteratively in cylindrical geometry using 

a point successive overrelaxation technique. Convergence between 

-6 

successive iterations was less than 10 after fifty iterations. 

The program was tested using ANL three group data. Flux shapes 
and energy spectra were determined for a typical fast reactor core 
and for a solid iron cylinder with a source at its center. The program 
was also used to determine criticality. 

Computation times were from one to ten minutes with less than 
15 OK words of core storage using the IBM 360/67. 
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I. INTRODUCTION 



A. BACKGROUND 

It has been found that the energy-dependent diffusion equation 
provides a fairly accurate description of the neutron behavior in fast 
reactor systems [1J. 

One way of conveniently using this equation is the now familiar 
multigroup method which was developed for reactor calculations by 
F. L. Friedman, A. M. Weinberg, and J. A. Wheeler [ 2 ] about 1947. 
Since that time, many solutions to the age-diffusion equation have 
resulted from the use of the multigroup method. Some analytical 
and semianalytical solutions have been carried out for a few group 
equations [3] but most of the work done has been numerical due mainly 
to the need for large numbers of energy groups, spatial points, and 
spatial regions . 

The numerical solutions to unique problems and classes of 
problems have evolved in the form of computer programs or codes . 

The degree of complexity and sophistication of these programs has 
roughly paralleled that of the digital computer. At this point in time, 
the literature abounds with numerous and varied codes and references 
thereto. Sangren [4] states that by 1959 there existed in. the United 
States more than 300 nuclear reactor codes for digital computers. 
Various groups have classified and categorized these codes [4] . In 
19 61 the American Nuclear Society and the Argonne National Laboratory 
established a code center for the purpose of collecting and disseminat- 
ing the multitude of codes that existed and for the future codes which 
were sure to be forthcoming [5], Around 1965 the European counter- 
part to this center was established £5] . 

Most of the more recent codes (those published after 19 65) 
examined in connection with this study are somewhat specialized 
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along the following lines: a) number of energy (lethargy) groups; 
b) number of spatial dimensions; c) number of regions; d) cross-section 
data input; e) geometry; f) number of mesh points; g) method of equation 
solution; h) finite difference approximation and i) many other features 
related to poisons, burn-up, perturbations, computer system, and others. 

B. OBJECTIVES 

This study was undertaken to provide a relatively simple-to-apply 
multigroup diffusion program for application to fast reactor cores and 
materials. In the future students at the Naval Postgraduate School 
will be able to use this program in conjunction with experimental 
neutron diffusion studies. In view of the large number of complex 
codes that exist and the difficulties related to modifying these codes 
for use on a specific problem and a specific computer system, it is 
anticipated that valuable research time may be saved by use of this 
program . 

Once the program was developed, initial testing would be done 
using a subcritical assembly. The first phase of this test would 
involve the solution of the three group diffusion equation in a typical 
fast core composition. From this, the energy and spatial dependence 
of the flux could be determined and comparisons made with the solution 
to the diffusion problem wherein space and energy dependence are 
treated as separable variables. 

The second phase of the initial test would involve the problem 
of the typical core assembly with a Pu-Be neutron source located at 
its center. This would then be modified to a solid iron cylinder with 
a centrally located Pu-Be neutron source. The flux shapes and 
energy spectra resulting from this study could then be used for 
comparison with experimental data resulting from local research. 
Experiments in this area are currently underway with applications 
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to fast reactor cross section data evaluation, materials and shielding 
reference £6] . 

Using the program to determine a critical composition and size for 
the typical fast reactor was another desired goal. In this area, a 
study was to be made on the effect mesh spacing has on criticality 
if any. 

Other objectives were to compare the energy spectra of the three 
and the sixteen group solutions and to ascertain any difference in 
the flux shapes and spectra by using the ANL sixteen group data 
and the LASL sixteen group data. 
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II. THEORY 



A. THE BOLTZMANN TRANSPORT EQUATION 

The conservation or balance of neutrons in a nuclear reactor can 
be completely described only by simultaneously specifying the distri- 
bution in space, time, energy and direction of motion of the neutrons. 
The general equation governing this balance is called the Boltzmann 
transport equation and may be written as 



PRODUCTION - LEAKAGE - ABSORPTION = 

a 

where -- is the time rate of change of neutron density. 
The neutron density is in general described as 



number of neutrons at time t in volume 
dr about r, whose speeds lie in dv about 
v and whose directions of motion (velocity 
vectors) lie in the differential solid angle 
angle dri about direction Jo. . 



and neutron flux is defined as 

fiCr, nr, _fL, t ) = Arn (r, nr, _n_,t ) 



ii- 1 



o 



II— 3 



Meghreblian and Holmes 1 describe Equation II- 1 verbally as 

number of neutrons at speed in dv about 
I v moving in direction -O. which appear 
per unit time from source in dr 



+11 



number of neutrons of direction -7L gained 
per unit time in dr from scattering 
collisions which scatter the neutron from 
all directions -ft - ' to -O- and all speeds 
to speed in dv about v 
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net number of neutrons moving in direction 
-FL at speed in dv about v that are lost 
per unit time by leakage through the 
boundaries of dr 



number of neutrons at speed in dv about 
moving in direction TL which are removed 
per unit time from dr by absorption 
collisions 



number of neutrons at speed in dv about 
moving in direction JrL which are removed 
per unit time from dr by being scattered 
into a new direction 

1 



=VI 



change per unit time in net number of 
neutrons of direction -Ti- in dr 



II— 4 



The first term in Equation II— 4 accounts for the rate at which 
neutrons with speed in dv about v and with direction in d-fi- appear 
from sources in dr . This term is then written as 



5 C Ej Ar, ) drdnrd_n_ . 



II— 5 



The second term of Equation II— 4 gives the total number of neutrons 
which appear in the differential element drdvd-pL per unit time 
with speed in dv about v and with direction dfn- about -Tu by 
being scattered from all possible speeds v and direction of motion -T i_ / 
The number of neutrons in drdvdxL which experience this scattering 
is 



22 s ( rv-') (j (r ; nrl st.\ "t ) dr drv'd sl' . II— 6 
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The probability that the scattering collision results in the neutron 
having speed in dv about v and direction a xl about -Tl is 
defined by 



5 C nr, xl j (V' sl' ) = 

probability that a neutron with an 
initial speed v and direction of 
motion xl , when scattered, emerges 
from the collision with a speed in 
dv about v and direction of motion 
in d XL about xi- . 



If Equations II— 6 and II— 7 are multiplied together and integrated over 
all possible v' and XL ! , the result is the total gain of neutrons in 

the differential element drdvdxL due to the scattering in of neutrons 
from higher speeds and different directions of motion written 



The third term of Equation II— 4 represents the net loss of neutrons from 
drdxL due to leakage through the elemental boundaries of dr and is 
equal to 



A rigorous derivation of Equation II— 9 is given in reference [l J . 

Term number four of Equation II— 4 is the absorption term and accounts 
for all neutrons lost from drdvdXL per unit time by absorption; it 
is given by 




dr dar'dx.' dXLdAk . 



II— 8 




II— 9 



( />Cr ) nr ) xv,t)drd nrdx\. 



II- 10 
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The rate at which neutrons are scattered out of drdvd-fL into new 
directions jrL' is the fifth term of Equation II— 4 and is written as 



T s Cat-) (pQr } nr, ) dr dru'dxi. . 11-11 

The final term of Equation II— 4 is the time-rate of change of the 
neutron density in drdvd Si. . This term can be written as 

§£ n(r> nr, _n_ , 1 : ) dr d r& dsL = 

(W ^ C. f j Ah, , t d f d fu' d _o_ 



By placing Equations II— 5 through 11-12 into the verbal Equation II— 4 
and dividing out the differential drdvd-TL the result is the Boltzmann 
transport equation for the distribution of neutrons of all energies in 
space and time 

“aF ft ^Cr, + X t O) (pCr, nr, IL/t) 

+ _TL* V (f) Cr, nr, _n_,t) ~ S Cr, nr, _fL,-t ) II- 13 

+ Xh jji' X s O') 5 Car, xL ; nr' } AOd/iFd-n-'. 



The second term in Equation 11-13 is the total loss of neutrons due to 
scattering out to new spaces and directions plus absorption. The total 
cross section is used. 

Equation 11-13 is perfectly general with the following two 
exceptions: a) the medium is assumed to be isotropic and homogeneous 
and b) neutron slowing down is due only to elastic-scattering collisions 
with nuclei . 
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of the gth group is defined by the integral 



^ 3 Cr) 




Cr/E) 



E 



dE 



11-17 



where the upper limit on the integral Eg-| is the lower boundary 
of the energy interval and the lower limit on the integral Eg is the 
upper boundary of the energy integral. The energy-dependent flux at 
the point F is $Cr>£') 

Within each group the diffusion of neutrons is described by the 
average diffusion coefficient. For simple diffusion theory this is 
written as 



d 3 = W — n - 18 

0 d ^-tr, 3 

where _ is the group macroscopic transport cross section and 

may be obtained from 






tp 







3 



11-19 



where the summation is over all elemental components and the group 

■7 

microscopic cross section cr_ in the absence of pronounced 

vc, g 

resonance effects due to capture and/or. scattering with the energy 
group is given by 




„Eg H 

E , u a- t z r , 9 CE) 0CE) dE 

P 3 ” ?SfE)dE 

J E 3 l 



11-20 
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B. THE AGE-DIFFUSION APPROXIMATION 

The age-diffusion approximation to Equation 11-13 is developed 
in reference [l] by expanding the Boltzmann transport equation in 
spherical harmonics. 

In the steady state the resulting first order approximation to 
Equation II- 13 is 



-DCu) V z 0Cr,u) + cu) 0 Cf,U) = SC f/u) 



11-14 



Since the change in the slowing down density of neutrons within the lethargy 
interval du must be balanced by the flow of neutrons scattered into 
du from elsewhere, the last term of Equation 11-14 may be written as 



r 



au 



2_ s u) 0Oi')chi' 



■' r T 1 f 

11- ic 



And Equation II- 14 may be rewritten substituting energy for lethargy 
variable as 



-D(E) ^ 2 0Cr ) E:) -v I a (E) 0 (r > E. ) = 

11-16 

SCr^ - ) 4r Z s C E.'-* E.") 0<e/) d.E.' 



Complete details of this analysis may be found in reference m- 

C. MULTIGROUP FORMULATION 

In forming the multigroup diffusion equation, the initial step is 
to divide the complete neutron energy (lethargy) range into N groups 
which may or may not be of the same width. The neutron flux 
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Neutrons may be removed from a particular group g through an 
absorption interaction described by the group absorption cross section 






°*> < 9 







eg 



+ in > 9 + ^ 



erjg 
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where 



g 

ln >3 



er,g 



= fission cross section of group g 

= capture cross section of group g 

= inelastic removal cross 
*5*9 section from group g 



=E 

K^3 



= elastic removal cross 

er>g-*K 

section from group g 
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The group macroscopic cross sections defined in Equation 11-21 above 
are obtained fiom the corresponding group microscopic cross sections 
for each distinct element or isotope in the reactor such that 



= E N Z cr^g , 11-23 

where q is the particular reaction and N equals the number of atoms 
of material z per cubic centimeter. 

The multigroup microscopic cross sections for each material are 
defined in reference [7] as follows: 




GROUP CAPTURE 



r 

Ccr^ (y (E)+ CE) + --E 0 CE)dE 



J £9 V(E)ciE 



11-24 
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f 5 9 



GROUP FISSION 



r-EgH 



e 9 



L cr n x >f CE) cj) CE) d£ 



p 3H (pC E ) d E 

Je 9 u 
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or 






GROUP SCATTERING (IN OR OUT) 

fEaH pt gH 

L. . J <?h, n CE-EO0fE) dE.de' 

Je =Ec^LO £ = Ecju j 



•K 



r 



E 9 H 






^(E)dE 
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where E T and E _ are the upper and lower limits of group g respectively 
gH gL 

and 0 (E) is the flux per unit energy interval at E. 

Neutrons may be added to a particular group g by a fission source 



Sf >9 ~ ^3 ^ ^ ^3 03 > 11 27 

3 " 1 

where % g is the fraction of fission neutrons born into group g such that 

r 

%o = % CE) dE , 

J J 631- 

<^^3 = JC N z Oo<r f ) 3 z 

z. 



11-28 

11-29 



and 



(Doj^ = 



, e 3 h 



e 3 l 



t) Z ( E-) o-p, f ( E") 0C&)dE. 



"E 3 h 

Egu 



0CE)dE 



11-30 
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The steady state diffusion equation for the gth energy group, where 



g represents any group between one and N, can now be written 




^3 0h (r) 

U - i 



Ds V 2 03 (r ) - r a , 3 0 3 (r) 

[9 — h)]0 9 (r)+ £ <J 't.Ch-^0 h (n 

C 0h (r) h = SgCr) 



11-31 



h- \ 



In Equation 11-31 the inelastic removal cross section from group g 
appears as the third term of the equation rather than as an implied 
part of the second term. This equation now represents a set of N 
coupled partial differential equations independent of time and applicable 
to one region of a reactor. The right hand side of the equation 
represents the number of neutrons at energies within group g coming 
from any extraneous neutron sources. 

D. SPACE-INDEPENDENT FORM OF MULTIGROUP EQUATION 

For this case, it is assumed that the extrapolation distance at the 
surface of the reactor is identical for all energy groups and thus is 
independent of energy. If this is true then all energy groups have the 
same spatial dependence. The flux is then written as 




11-32 



where g is the magnitude of the flux in the gth energy group, and 
the function F(r) satisfies the equation 



V 2 F (r ) + B 2 F(r ) = 0 



11-33 



2 

the factor B is the square of the first eigenvalue of Equation 11-33 
and is called the buckling. 
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By substituting Equations 11-32 and 11-33 into Equation 11-31 
expanding the resulting equation in a series of eigenfunctions, and 
applying the orthogonality condition, the result will be a set of N 

linear algebraic equations in N unknowns written 

* 

-[D 3 B 2 + Xa l3 +13 X Cg-M ] 

9-» X h >9 N ^ J 

+ 5zf h + X 3 = So . 

h=i h=i 



This set of equations is easily solved by Cramer's Rule. 

E. EXPLICIT SPATIAL DEPENDENCE 

In order to solve the set of Equations 11-31 by numerical methods, 
the Laplacian operator in the first term must be replaced by its equiva- 
lent finite difference approximation. The resulting equation for a 
typical energy group g in cylindrical geometry is 

Dg [-p { ( r-t-i , z, ©) ~ 20 3 (r,z, 8 )+ 0 g(.r-i,z,e)} 

+ 2Fh£ 03( r +l,z,e)-0 3 (r-i ) z,6)3 +T? { 0 g Cr,z-i,e) 

(r,z,©) +09(r J z+i > e)3 + - f 2 j i £ 0 3 (r,z, e + 0 

-209(r 1 z,e)+09(r J z ) e-i)3] 'Za.j 0 3 Cr,'z,e) 

-[|lZC9-D] 9Vr,z,e>+ S Zth-^pfhCnz.e") 11-35 
+ Xg ST'Oh'Ej.v, 0h(r,z,e') = S 9 o-.z, s') . 

h-l 



Appendix A contains the details of the derivation of the finite difference 
approximation and the definition of the symbols used in Equation 11-35. 
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Throughout this investigation, circular summetry was imposed. 
Thus the third term in Equation 11-35 is zero. For mesh points not 
on the longitudinal axis (See Appendix A for finite difference operator 
used on this axis.), Equation 11-35 becomes 

Dg 0 3 Cr + i , - 2 .) - 2 0s (r,z2) + 03 Cr-i,zA 1 
+ 27h l 09 Cr + 1 > + { ^gCr, z+O 

- 2 ^ 3 (r ; z) + 09 C r > z.-\) ^ - Zcx,g 0 g (r ,-z*) n . 

N <3-> 

-[ |2 0g cr,z) + 22 Zct'-<a')0h(r ) z-) 

h>9 N IN--I 

+ X 3 22 51 f)K ph(r,z) = S 3 (r > x.') . 

h-l 



F. MATRIX FORM OF MULTIGROUP EQUATIONS 

For the purpose of numerical analysis and computational ease. 
Equation 11-36 is reduced to the matrix form. In Appendix B, this 
equation is first expanded to fit the available sixteen group data, 
then reduced to a form more suitable for matrix notation. Finally, 
it is reduced to the form 

A l 0g (r-i, -z.)+ A 2 0g (r,z) + A3 0g Or 4- \,'z') 
+A40 9 (r, z-0 +A5 0g(r, z+i) = 

B-4 

4 - AS, 0, Or/ 2 ) 4 - • • • 4- A S 9-1 09-1 C r,z) 

4- AFg 0g (r,z') 4- • • • 4- AF l8 0i6 Cr , z:2) 
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To reduce this set of multigroup equations to matrix form, the 
spatial dependence is considered first. 

The flux vector of group g is defined as 



0 ) 1 ) 







09 ( 1 ) tv\ ) 



11-37 



where ^>gCL ; 2- + M) and 0g(2~*L ) M) are on the extrapolated 
boundary and are set equal to zero. The fluxes 0CI>X) and 0(X)O 
are set equal to the fluxes 0C3.X) and 0CX, 3) respectively 
in order to satisfy the normal derivative boundary condition. See 
Appendix C for mesh point numbering scheme. 

The matrix of coefficients Ag is written 



O AA O' 
O O AA 



A s - 



O 



O AIA2A30--0 OA50--0--0 
O OA1A2A3--0 OOA5O--0 



O O 



A4 0--AI 



A5 



) 11-38 



The terms on the right hand side of equation B-4 are called the 
source terms and Wg is the source vector for group g. It con be 
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represented in matrix form by 



3 -' '8 

Wg = X - 0h + AF h 0b + Sg 0b 

Vr=l h-g+> 



where the flux vector is given as 




0b 0,1) 

0h (2,\) 

0h(L, 0 
0 h (L ; M) 



the scattering matrix ASh 



O • • • O O ASh O • • O o O 

O • - • O O O ASv, -O Q O 



AS 



h 



O • • • O O o O • * • ASv> 

the fission matrix AFh 



O - • • O OAF h O - ■ O O O" 

O • • • O O o AFh - -O O O 



AFh - 



0---0 O o 



O- • • -AF-, 



11-39 



11-40 



11-41 



11-42 
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and the extraneous source matrix Sg 



O * * * * O O Sg o • • • o o o 

O • • • • O O O • • O O o 




o . - . . o o 

L 




11-43 



The result is then a set of equations, one for each spatial mesh 
point for the energy group g. These may be written as 






11-44 



Next the equations are cast according to the group dependence. 
Define czS as 

4 




11-45 



A as 



A, O ' • 

O A z • 



O • 




11-46 
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and W as 



w 



w, 



w N 



The entire set of equations then takes the form 



11-47 



A (f> = W 



11-48 



G. NUMERICAL SOLUTION OF THE MULTIGROUP EQUATIONS 

The basic equation 11-48 can be solved by various direct or 
iterative methods. Direct methods generally are based on either 
Gaussian elimination with pivoting, or on triangular decomposition of 
the matrix of coefficients. Iterative methods include the Jacobi 
method, unextrapolated Liebmann or Gauss-Seidel method, various 
forms of the extrapolated Liebmann or successive orerrelaxation 
method and others. 

In this study, up to 25,000 algebraic equations were expected as 
a result of a sixteen group problem in two dimensions with 1600 mesh 
points. For this large number of equations, iterative methods are 
generally used rather than direct methods £8, 9]. Smith [9] states 
that one reason for this is that the iterative methods are more efficient 
in that they take advantage of the sparse matrix of coefficients with 
its numerous zero elements . 

In this study, the point successive overrelaxation method was 
used (See Appendix D) . 
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First the reactor materials and energy group structure are chosen. 

Next the cross section data is computed or as was the case in this 

study, taken from some reliable source [7 ] and used to compute the 

appropriate coefficients. Then initial estimates of the group fluxes 

designated 0g° at all mesh points in the reactor except those lying 

on the extrapolated boundaries are made. These values are then used 

along with their respective coefficients to compute the initial source 
© 

vectors Wg . The flux in group one at all mesh points in the 
reactor is then corrected by solving the set of equations 

A, 0,° = W,° 11-49 



Using the point successive overrelaxation method, this is done 
by solving the first of Equations 11-49 for the central mesh point, 
i.e., 0,C2,2) and then solving the remaining equations point- by- 

point, row-by-row, at all other points no on the extrapolated boundary. 
The solution is denoted 0, . The remaining group equations are then 

similarly solved yielding the vector 0 ' where 





■) 




11-50 



which represents the first iteration. 

This vector is then rescaled such that 
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which indicates that the flux vector on the zeroth iteration integrated 
over all spatial dimensions and all energy groups is equal to the flux 
vector after the first iteration integrated over all spatial dimensions 
and all energy groups multiplied by a scale factor [10] . When the 
process is carried out, it can be seen that after a few iterations the 
scale factor will approach a constant value. The reason for this 
as described by Lamarsh 1 11] is that each iteration is actually 
equivalent to one cycle of the chain reaction. The zeroth iteration 
of the fluxes 9 then can be thought of as a set of initial conditions 
or the flux state at time t = 0. Eventually the flux will approach the 
funadmental eigenfunction as the higher harmonics die out. The 
fundamental will then vary with the number of iterations (time) 
increasing or' decreasing depending on whether the system is 
subcritical or critical. The quantity $ can clearly be considered 
a measure of the criticality of the system. If S is less than unity 
then the group fluxes throughout the system are decreasing with each 
iteration and the system is subcritical and if is greater than 

unity then the fluxes are increasing and the system is supercritical. 

"S is in fact the inverse of the effective multiplication factor K 



The iteration process described above is referred to as the inner 
iteration cr flux iteration. If a critical assembly is desired, then 
an outer iteration must be performed wherein adjustments are made to 
the composition or size of the reactor in order to make the effective 
multiplication factor converge to unity. 

Convergence of the flux vectors on the inner or flux iterations 
is determined by taking the absolute value of the difference between 
the flux at all mesh points and in all energy groups and comparing it 
with some desired value £ [ , such that 



eff 
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Convergence of the outer or criticality iteration is determined by 
subtracting the scale factor "S from unity so that 



1.0 - 'b 
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where €2. is an arbitrary convergence criterion. 
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III. DESCRIPTION OF COMPUTER PROGRAM 



A. INPUT DATA 

In this section the input data is read in via punched cards as 
illustrated in Appendix E and printed in the output for verification. 

The input data along with their symbolic program designators include 
the following: the energy dependent source SS (I); the maximum values 
of the initial flux estimate S (I); the energy intervals DELE(I); the number 
of atoms per cubic centimeter H(I) of each element or isotope; the 
microscopic cross sections for each element or isotope for inelastic 
scattering SIGIN(I,J,K) , transport SIGTR(I,J), fission SIGF(I,J), 
elastic removal SIGER(I / J) / and capture SIGC(IJ); the normalized 
neutron fission spectrum CHI(I); and the average number of neutrons per 
fission NU(I,J). 

B. OUTER ITERATION 

If a critical composition is the desired result, the program returns 
to this section after the iteration has converged (See SAMPLE COMPUTER 
PROGRAM page 119). The iteration on composition then begins by 
adding to the volume of the reactor an incremental volume of uranium 
235. An equal volume of uranium 238 is subtracted and new values for 
N^S H(l) and N^ H(2) are computed. On the zeroth outer iteration 
this section of the program is by-passed. 

C. COMPUTATION OF MATRIX OF COEFFICIENTS 

Input data is used in this section to compute the elements of the 
coefficient matrix, the scattering matrix, and the fission source matrix. 
In the first part of this section, the various microscopic cross sections 
for each reactor component are multiplied by their respective atomic 
abundances N . The values for each particular reaction (for example, 
capture or inelastic scatter from g to g + 1) are then summed. This 
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is done for each energy group. The summations are then combined 
with other input data and the coefficients are computed for each 
group equation. 

The equation denoted by statement number 17 in the sample 
program computes the diffusion coefficients for each group based on 
the simple diffusion theory approximation. The absorption coefficients 
for each group are computed in the equation denoted by statement 
number 175 by summing the contributions from all the absorption 
or loss reactions. The equation denoted by statement number 18 
computes the fission and down- scattering coefficients from all 
higher energy groups. This is accomplished by starting at the 
highest energy group in which the reactions originate and computing 
its contribution to all lower energy groups at the end of the slowing- 
down process. The fission source coefficients for a particular group 
and all lower energy groups are computed by the equation given in 
statement number 21. The above coefficients are completely defined in 
Section II and Appendix B. 

D. INITIAL FLUX ESTIMATE AND BOUNDARY CONDITIONS 

In the cylindrical geometry the flux shape is initially estimated 
as the product of the relative group magnitude factor S(I), a zero 
order Bessel function, and a cosine function. The factors S(I) were 
determined for this study as suggested by Clark and Hansen [10] . 

The flux on the exterior boundaries is set equal to zero. At the 
interior boundaries the flux at the first mesh points on either side of 
the longitudinal and mid-plane axes are set equal, indicating that the 
flux has its maximum value along these axes. At the end of this section 
the initial flux estimate is printed out for verification and comparison 
with the final flux values . 

The program is designed to by-pass this section on all outer 
iterations . 
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E. INITIAL FISSION SOURCE INTEGRATION 

The initial flux estimate is used in this section to compute the 
fission source vector for the zeroth inner iteration. The program 
performs the following triple integration: 



r E 



O vJ< 



{[ 0<E.,r,z)] %(E) (E)jdEdrd 



z 



hi- 1 



by dividing the energy spectrum into discrete intervals DELE(I) 
corresponding to the energy groups (for this study the upper boundary 
of the highest energy level was taken as 10 Mev for use in this calcu- 
lation) and dividing the area into finite intervals corresponding to the 
mesh spacing GS. A summation process is then carried out. 

F. COMPUTATION OF SOURCE VECTORS 

The components of the vector W for each inner iteration are computed 
by summing the products of the fission source or the combined fission 
source plus down-scattering source coefficients and the scaled flux 
values from the previous iteration. On the zeroth iteration the initial 
flux estimate is used with a scale factor of unity. The resulting numbers 
SUMR(I,R,Z) and SUML(I,R,Z) are the space dependent components of 
each group source vector Wg. 

G. SOLUTION OF THE DIFFUSION EQUATION 

This section consists of one outer energy group loop with several 
space loops within it. Each group equation beginning with the highest 
group is solved at all spatial points starting at the longitudinal axis r = 0 
and proceeding column by column to within one mesh point of the extra- 
polated boundaries r = R-l and z = Z-l. 
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The first space loop will solve the diffusion equation at the 
geometric center of the reactor, i.e. , point (2,2) (See Appendix C, 
Figure C-l.) if an extraneous point source is located in this position. 
This particular configuration was used in this study; however, extrane- 
ous sources with other spatial distributions could also be used with 
the appropriate changes to this section of the program. 

The second space loop involves the solution of the diffusion 
equation at all points on the longitudinal axis r = 0 but not occupied 
by an extraneous source. A special finite difference approximation 
for points on this axis is necessary due to a singularity caused by the 
radius being equal to zero (See Appendix A, page 66). 

The last spatial loop solves the diffusion equation at all remaining 
mesh points with the exception of those located on the exterior 
boundaries . 

Within each of these spatial loops, the flux value at each mesh 
point from the previous iteration is unsealed (TPHI1) and compared 
with the flux value at each mesh point in the present iteration to 
determine convergence as defined in Section II , page 27. The 
largest absolute value of these differences, 



TE.ST2 = DABS ( TPH1Z - TPHlO , III- 2 

is stored as the quantity SAVE and may be printed with the output as 
desired. The scaled flux values of the previous iteration SPHI1, 

SPHI2, SPHI3 are used in the successive point overrelaxation scheme 
for solving the group equations at each mesh point. 

H. NTH FISSION SOURCE INTEGRATION (UNSCALED FLUX) 

In this section the triple integration process, Equation III— 1 , is 
performed with the new but as yet unsealed flux values. 
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I. COMPUTE SCALE FACTOR AND RENORMALIZE FLUX 

Once the fission source integration has been computed for the 
N-l iteration FLUXI(NN-l) using the scaled or renormalized flux 
values and for the Nth iteration FLUXI(NN) using unsealed flux values, 
the scale may be determined from 

SCALF = DSQRT(FLUXI(NN-1)/FLUXI(NN)) III- 2 

The flux values at all mesh points and in all energy groups are then 
renormalized by multiplying by this scale factor. The scale factor 
SCALF is the principal eigenvalue of the system and its reciprocal is 
the effective multiplication constant SKEFF as explained in Section II, 
page 



J. NTH FISSION SOURCE INTEGRATION (SCALED FLUX) 

At this point the integration process of Section III (E) is performed 
using the renormalized flux. On the subsequent iteration this will 
become the quantity FLUXI(NN-l) in the numerator of Equation III— 2 . 

K. END OF INNER ITERATION 

This statement (1801) is the controlling statement for the inner 
iteration. It can make the program exit to either the output section 
or the outer iteration loop depending upon whether or not a critical 
reactor is desired. The controlling factor may be either a specific 
number of inner iterations if experience has shown that the convergence 
criterion as defined in Section II, page 27 will be satisfied or another 
possibility would be to stop the inner iteration once the convergence 
criterion is met regardless of the number of inner iterations. The 
second method might be preferable once it had been determined that 
convergence to a reasonable criterion, i.e. , 6, would not take an 
excessive amount of computer time. 
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L. OUTPUT SECTION 

The output is very flexible and may be easily changed to fit a 
particular problem or set of problems. Typical output data includes the 
following: scale factor, effective multiplication constant, flux 

values at various mesh points and various inner or outer iterations, 
convergence criterion SAVE ( € t ) and overrelaxation factor w. In 
addition, graphical output in the form of graphs of flux shape compared 
with the cosine function and zero order Bessel function at various 
locations may be plotted and the convergence criterion SAVE versus 
iterations. Many other possibilities exist. 

M. START NTH OUTER ITERATION AND TEST FOR CRITICALITY 

In this section the scale factor is checked out for criticality 
by comparison with an arbitrarily selected value ( .0005 was chosen 
for this study.) such that 

i 1 . 0 — SCALF | £ .0005 III-3 

If the above condition is not satisfied, then another outer iteration 
is made. The number of outer iterations may also be arbitrarily 
limited in this section. Once this limit, for example, ten outer 
iterations, is reached, the program will stop even though criticality 
has not been reached. 
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IV. RESULTS AND PROGRAM VERIFICATION 



If in the paragraphs that follow modifications to the basic core 
described below have been made, these modifications will be so stated. 
The results described and tabulated in this section were based on a 
typical homogeneous fast reactor core composed of the following 
materials with their respective volume fractions: 



U 

U 



235 

238 



Na - 



Fe 

Pu 



239 



2 % 

30% 

33 % 

27% 

8 % 



This composition was chosen to coincide approximately with the 
RAPSODIE reactor core in reference [ 12] . The core dimensions were 
36 centimeters in radius and 72 centimeters in height. Mesh spacing 
was two centimeters in both the radial and axial directions . Three 
energy groups were used with the cross section data taken from 
ANL 5800 [7]. 



A. DETERMINATION OF OPTIMUM RELAXATION FACTOR 

The experimental method used for determining w ^ consisted of 
two basic steps. First, an analytical estimation was obtained based 
on formulation developed by Frankel [13] and Starling £14] . 

Frankel [13] derived the following equations based on the 
solution to Laplace's equation in a bounded region R in rectangular 
coordinates 



w opt = 4 OC 



IV- 1 



f 
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♦; 




where ot is the smaller root of 



Or 2 t 2 - A ol ■+■ I - O 



IV- 2 



and t is given by 

"t = cos + COS IL IV- 3 

where p and q are the number of mesh points in the x and y directions. 
Smith [9j used this same method to solve the similar Poisson's equa- 
tion in rectangular coordinates in a bounded region R. Starling [14] 
has developed an equation for t applicable to Laplace's equation in 
cylindrical coordinates and given by 

7 - ~ COS ^ + To C ^ <5 IV- 4 

where Z and R are the height and radius of the cylinder respectively 
and h is the mesh spacing. 

The table below gives the results found for these analytical 

solutions for w using a symmetric cylindrical section 36 centimeters 
opt 

in radius and 36 centimeters in half-height 



METHOD TOR FINDING "t" 


w . 
opt 


REMARKS 


Eq. IV- 4 


1.7888 


Z = 72 cm. (full height) 


Eq. IV- 4 


1.7259 


Z = 36 cm. (half height) 


Eq. IV- 3 


1.7294 


p = q = 20 (no. mesh pts . + 1) 


Eq. IV- 3 


1.704 


p = q = 18 (no. mesh pts.) 



The second step in the experimental process was to try these values 

in the actual program and graph the convergence criterion versus the 

number of inner iterations for these values of w along with other 

opt 

reasonable estimates to determine the best one. Figure IV- 1 shows the 



37 




N0183XI ao 30N3SH3AN03 



38 



ITERATIONS 



data which resulted from this study. It can be seen that for the values 



of w ^ equal to 1.72, 1.73 and 1.74, a convergence criterion value 
opt -6 

of less than 10 is satisfied within fifty iterations. Based on this 

data , it was concluded that Equation IV- 4 with Z equal to one-half 

the height of the cylinder, due to symmetry, gave the best value of 

w for this program. The program was modified to compute w 
opt opt 

automatically based on Equation IV- 4. 

The above results were reasonable in view of the fact that 
Equation IV- 4 was developed for use with Laplace's equation in 
cylindrical coordinates . The multigroup neutron diffusion equation 
used in this study in cylindrical geometry is of a similar form. One 
of the major differences is that the dependent variable 0 is a 
function of both space and energy. 



B. FLUX SHAPES FOR A TYPICAL NONCRITICAL CORE 

For the bare reactor described at the beginning of this section, the 
initial flux shape was estimated to be the product of the relative group 
weighting factor S(I), a cosine function and a zero order ordinary 
Bessel function of the first kind. This corresponds to the solution of 
Laplace's equation in a finite cylinder. Figure IV- 2 shows the final 
flux shape in the axial direction at the center of the reactor after 
260 inner iterations with a convergence criterion of less than 10 ^ 
compared with the cosine shape at the same location at zero iterations . 
The magnitudes were normalized to unity so that the shapes could be 
easily compared. It was found that the deviations were in the order of 
1 to 15 per cent with the per cent deviation increasing from the center 
outward. All three energy groups showed exactly the same normalized 
shapes after an equal number of iterations. The magnitudes of the 
group fluxes are unknown since they depend on the power level at 
which the reactor is being operated. 
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NORMALIZED FLUX 



AXIAL FLUX SHAPE 



IN A TYPICAL SUBCRITICAL CORE 




Figure IV- 2 
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The same comparison has been made in Figure IV- 3 between the 
flux in the radial direction and the assumed initial Bessel function. 

A very slight deviation in shapes was noted here also indicating the 
close similarity between the diffusion equation and Laplace's 
equation . 

Examination of the data revealed that the final flux shapes were 

in effect established completely to seven decimal places after some 60 

inner iterations at which time the convergence criterion SAVE was less 
-7 

than 10 

The program was run using a constant as the initial flux shape 

estimate rather than the cosine and Bessel function. The resulting 

final shapes for all three groups were exactly the same as those in 

Figures IV-2 and IV- 3 after 90 inner iterations. Convergence was 

somewhat slower using this initial estimate. For example, after 50 

-5 -7 

iterations, SAVE was equal to 3 . 2 X 10 versus 3.7 X 10 for an 
initial estimate based on the cosine, Bessel function product. 

C. PROGRAM TEST WITH EXTRANEOUS POINT SOURCE 

The program was run with a simulated extraneous Pu-Be point 
neutron source at the geometric center of the reactor. Two core 
compositions were used in these tests. 

The first core was made up of the typical composition described 
at the beginning of this section. The flux shapes in the radial and 
axial directions for the three groups were plotted on the same graph 
with the Bessel functions and cosine respectively. Figures IV- 4 and 
IV- 5 show the results of these comparisons after 90 inner iterations. 
Group one showed the largest deviation with the most pronounced 
effect nearest the source. This result was due to the normalized source 
neutron contribution to group one. Eighty- seven per cent of the neutrons 
from the Pu-Be source are born with energies lying in group one 
(1.35 Mev —►co ) as defined in ANL 5800 [ 7 D . 
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NORMALIZED FLUX 



RADIAL FLUX SHAPE 



IN A TYPICAL SUBCRITICAL CORE 




Figure IV- 3 
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NORMALIZED flux 



RADIAL FLUX SHAPES 



IN A TYPICAL SUBCRITICAL CORE 
WITH A CENTRAL POINT SOURCE 




Figure IV- 4 
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AXIAL FLUX SHAPES 



IN A TYPICAL SUBCRITICAL CORE 
WITH A CENTRAL POINT SOURCE 
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The second core was a solid iron cylinder with the same 
dimensions as the typical fast reactor core. The Pu-Be source 
was located at the center of the cylinder. The radial flux shapes for 
groups one and two were plotted in Figure IV- 6 and the axial flux 
shapes for these two groups were plotted in Figure IV- 7. It was noted 
that for these two groups the flux decayed very rapidly as the distance 
from the source increased; therefore, the ordinates of Figures IV- 6 
and IV- 7 are logarithmic. Figures IV- 8 and IV- 9 show the radial and 
axial flux for group three. In these graphs the flux shapes are com- 
pared with the cosine and Bessel functions after 90 inner iterations. 

The three group flux spectrum at the center of the two assemblies 
are shown in Figure IV- 10. A relatively flat spectrum was observed for 
the core containing the fissile material while the solid iron core showed 
large relative differences among the three group fluxes. 

D . CRITICALITY TEST 

Test of the program up to this point included only the inner iteration. 

Criticality was achieved by two methods of outer iteration. 

The first method consisted of performing sufficient inner iterations 

for convergence. Seventy iterations were used resulting in an of 

-9 

less than 10 for a mesh size of two centimeters. An outer iteration 

235 

was then carried out wherein .5% by volume of U was added before 
each outer iteration as described in Section III-B, This cycle was then 
repeated until the criticality criterion ( | 1.0 - SKEFF l<.Ooo5 ) 

was satisfied. The optimum overrelaxation factor remains the same 
for each outer iteration. 

Four different mesh spacings were used in this section of the 
program analysis to determine the effect, if any, mesh size had on 
criticality and flux shapes. Table IV- 1 contains the data resulting 
from this study. The central flux shapes in the radial and axial 
directions at criticality for all mesh sizes tried were nearly coincident 
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NORMALIZED FLUX 



RADIAL FLUX SHAPE 
IN AN IRON CYLINDER 
WITH A CENTRAL POINT SOURCE 
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NORMALIZED FLUX 



AXIAL FLUX SHAPE 



IN AN IRON CYLINDER 
WITH A CENTRAL POINT SOURCE 




Figure IV- 7 
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normalized flux 



RADIAL FLUX SHAPE 
IN AN IRON CYLINDER 
WITH A CENTRAL POINT SOURCE 




Figure IV- 8 
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3)V 



NORMALIZED F LUX 



AXIAL FLUX SHAPE 



IN AN IRON CYLINDER 
WITH A CENTRAL POINT SOURCE 




HEIGHT from center ccm.) 



Figure IV- 9 
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FLUX 



FLUX SPECTRUM 
AT CORE CENTER 
WITH A CENTRAL POINT SOURCE 
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Figure IV- 10 



S-l 



10 ° 10 ' 



50 



MESH SPACING = .75 CM. 



VOL. PCT . 
Xj235 


K rr AFTER 70 INNER 
pff 

ITERATIONS 


OUTER 

ITERATION 


EPSILON 1 


2 


.99464 


0 


1.8 X 10“ 4 


2.5 


.99567 


1 


3.9 X 10" 5 


3 


.99668 


2 


3.7 X 10" 5 


3.5 


.99768 


3 


3.5 X 10" 5 


4 


.99865 


4 


3.3 X 10" 5 


4.5 


.99968 


5 


3.3 X IQ" 5 


CORE STORAGE 270K 


COMPUTER TIME 33:26 




MESH SPACING : 


= 1.5 CM. 




VOL. PCT. 


K AFTER 70 INNER 


OUTER 




U235 


6ff ITERATIONS 


ITERATION 


EPSILON 1 


2 


.98952 


0 


4.7 X 10" 8 


2.5 


.99156 


1 


1.7 X 10~ 8 


3 


.99355 


2 


2.5 X 10" 8 


3.5 


• .99550 


3 


3.4 X UT 8 


4 


.99741 


4 


4. 2 X 10 -8 


4.5 


.99913 


5 





CORE STORAGE 148K COMPUTER TIME 8:34 



MESH SPACING = 2.0 CM. 



VOL. PCT. 
U 2 3 5 


K AFTER 70 INNER 
eit ITERATIONS 


OUTER 

ITERATION 


EPSILON 1 


2 


.98631 


0 


1. 1 X 10" 9 


2.5 


.98898 


1 


6.0 X 10" 10 


3 


.99159 


2 


6.0 X 10“ 10 


3.5 


.99415 


3 


9.0 X 10" 10 


4 


.99665 


4 


1.2 X 10 -9 


4.5 


.99910 


5 





CORE STORAGE 126K COMPUTER TIME 5:40 



MESH SPACING =4.0 CM . 



VOL. PCT. 
Tj235 


K AFTER 70 INNER 
e “ ITERATIONS 


OUTER 

ITERATION 


EPSILON 1 


2 


.97542 


0 


lo-io 


2.5 


.98036 


1 


10 -10 


3 


.98517 


2 


io-io 


3.5 


.98985 


3 


10-10. 


4 


.99442 


4 


10“ 10 


4.5 


.99 


5 





CORE STORAGE 100K COMPUTER TIME 1:41 



Table IV- 1 
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with the Bessel function and cosine respectively. Table IV- 2 shows a 

maximum deviation of .71 per cent using a mesh size of two centimeters. 

235 

Criticality was reached at the same core composition, i.e. , 4.5% U , 
for mesh spacings of .75 centimeters, 1.5 centimeters and 2 centimeters. 
With a mesh size of four centimeters, criticality was not achieved 
at this composition. 

The second method for arriving at a critical assembly consisted 

of fixing the core composition and varying the core dimensions in the 

outer iteration. A new optimum overrelaxation factor was computed 

for each outer iteration since the geometric parameters that were used 

to compute w were changing. Three compositions corresponding 
opt 

to 18.9%, 20% and 27% enrichment were used. The results of this 
test are plotted in Figure IV- 11 and compared with the bare core criti- 
cality curve of the RAPSODIE reactor in reference [11]. Differences 
noted are believed to be due to the following: a) compositions are not 
exactly the same; b) different cross section data was probably used 
though this was speculative since reference [11] did not state the 
source of their cross section data; and c) different programs were used 
to solve the multigroup diffusion equations. PRODII program was used 
for the RAPSODIE reactor. The factor above having the most influence 
was probably the large volume of iron. The absorption cross section 
of iron is fairly large and as seen in Table IV- 3 is occupying a volume 
which in RAPSODIE contained 4% molybdenum, an unspecified quantity 
of iron in the form of stainless steel, and possibly some void space. 

The absorption cross section of these components is likely to be less 
than pure iron. 

The comparison of these curves shows that while a much larger 
critical reactor, assuming the same enrichment, would result from 
using the program developed in this study, the shapes of the curves 
are very similar. The fact that as enrichment was increased, the core 
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MESH POINT 
LOCATION 


INITIAL FLUX 
GROUPS 1, 2, 3 


FINAL FLUX 
GROUPS 1,2, 3 


PER CENT 
DEVIATION 










(2,2) 


1.000000 


1.000000 




(4,2) 


.982230 


.981311 


.090 


(6,2) 


• 929858 


.928161 


. 182 


(8,2) 


. 845676 


.843389 


. 270 


(10,2) 


.734131 


.731515 


.356 


(12,2) 


.601101 


.598447 


.441 


(14,2) 


.453505 


.451134 


.523 


(16,2) 


.298939 


.297153 


. 597 


(18,2) 


. 145212 


. 144270 


.648 


(20,2) 


• 0.0 


0.0 












(2,2) 


1.000000 


1.000000 


— -- 


(2,4) 


.984808 


.983927 


. 089 


(2,6) 


.939696 


.938015 


. 178 


(2,8) 


.866030 


.863707 


.268 


(2,10) 


.766040 


.763312 


.356 


(2,12) 


.642787 


.639923 


.445 


(2,14) 


.500000 


.497327 


.534 


(2,16) 


.342020 


.339888 


.623 


(2,18) 


. 173646 


. 172411 


.710 


(2,20) 


0.0 


o . o - 





PER CENT DEVIATION = 



INITIAL FLUX - FINAL FLUX 
INITIAL FLUX 



X 100 



Table IV- 2 



BARE REACTOR 
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RAPSODIE 


18.9% ENRICHED 


20% ENRICHED 


COMPUTER 


IBM 650 


IBM 360/67 


IBM 360/67 


COMPUTER CODE 


PROD II 


THIS STUDY 


THIS STUDY 


NUMBER 
OF GROUPS 


? 


3 


3 


CROSS SECTION 
DATA 


? 


ANL 3 GROUP 


ANL 3 GROUP 


COMPOSITION 
VOL. PERCENT 








URANIUM 


28% 


32% 


32% 


PLUTONIUM 


8% 


8% 


8% 


SODIUM 


33% 


33% 


33% 


MOLYBDENUM 


4% 


— 




STAINLESS STEEL 
IRON 


? 


27% 


27% 



Table IV- 3 



size and hence the critical mass was decreased indicates a correct 
trend. The flux shapes found using this outer iteration scheme were 
also nearly coincident with the cosine and Bessel function curves. 

Table IV- 4 shows an example of the per cent deviations experienced 
for a critical core assembly 32 centimeters in radius and 124 centimeters 
in height. 

E. TESTS WITH SIXTEEN GROUP DATA 

Attempts to modify the program to solve a sixteen group problem 
were not successful. Both under- and over-relaxation factors were 
tried but satisfactory convergence of the inner iteration could not 
be achieved after over thirty minutes of computer time . An under- 
relaxation factor of . 1 was used and after 310 iterations and 34 
minutes, the convergence criterion was only of the order 10 
Overrelaxing with a factor of 1.726 yielded a highly oscillatory flux 
behavior in all energy groups. All group fluxes were initially estimated 
to have maximum value of ten. The magnitude of oscillations recorded 
in the output varied between plus or minus eighty for the higher energy 
groups to plus or minus thirty for the lower energy groups with no 
noticeable damping after 310 iterations and some 34 minutes of 
computer time. No further tests were conducted. 
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MESH POINT 
LOCATION 


INITIAL FLUX 
GROUPS 1, 2, 3 


FINAL FLUX 
GROUPS 1,2,3 


PER CENT 
DEVIATION 










(2,2) 


1.000000 


1.000000 


— 


(4,2) 


.977535 


.977479 


.006 


(6,2) 


.911646 


.911554 


.010 


(8,2) 


.806757 


.806654 


.013 


(10,2) 


.669888 


.669786 


.015 


(12,2) 


.510080 


.510013 


.013 


(14,2) 


.337808 


.337782 


.008 


(16,2) 


. 164121 


. 164154 


.020 


(18,2) 


o 

♦ 

o 


o 

o 


• 


(n 


1.000000 


1. 000000 


— 


(2,4) 


.994868 


.994805 


.006 


(2,6) 


.979525 


.979403 


.012 


(2,8) 


.954141 


.953952 


.020 


(2,10) 


.918959 


.918717 


.026 


(2,12) 


.874343 


.874059 


.032 


(2,14) 


.820767 


.820438 


.040 


(2,16) 


.758757 


.758407 


.046 


(2,18) 


.688969 


.688603 


.053 


(2,20) 


.612101 


.611742 


.058 


(2,22) 


.528959 


.528615 


.065 


(2,24) 


.440393 


.440076 


.072 


(2,26) 


.347303 


.347033 


.078 


(2,28) 


.250656 


.250440 


..086 


(2,30) 


. 151424 


. 151289 


.089 


(2,32) 


.050646 


.050600 


.091 
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V. CONCLUSIONS 



This study was undertaken to develop an easy to apply computer 
code for numerically solving the multigroup diffusion equations. The 
program was written in the FORTRAN IV language and all work was done 
on an IBM 360/67 computer. Applicability to other computer systems, 
while not a consideration of this study, is assumed to be practical. 

The unique feature of this program is that it will be readily 
available to and easily used by students at the Naval Postgraduate 
School. It is hoped that this study will complement future local 
experimental studies in fast neutron diffusion and also enable the 
student researchers to conserve valuable time which would be 
expended attempting to verify their experimental data with a myriad 
of complicated and unfamiliar multigroup diffusion programs. 

The point successive overrelaxation method of solution was 
successful with the three group data on the basis of the convergence 
criterion chosen, i.e. , 

| 0 NM CE,r>z) ~ 0 N ( E,r, Z)| < IO~ 6 

all 
E- ,r, -z. 

Convergence of the three group criticality problem with five outer 
iterations took five minutes and forty seconds and 126K-words of 
core storage using a mesh size of two centimeters. Criticality was 
achieved with smaller mesh sizes, i.e. , .75 centimeters and 1.5 
centimeters, but computation times were considered excessive especially 
for .75 centimeter mesh size which took about 33 minutes of computer 
time and over twice as much core storage space. 

Tests showed that the code was easily adaptable to three group 
solutions for typical fast reactor cores and reactor material assemblies 
with extraneous point sources. Flux shapes and energy spectrums from 
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these tests should provide a sound basis for future experimental flux 
measuring studies in fast reactor materials. 

The close coincidence of the critical flux shapes with the cosine 
and Bessel functions indicates that for a three group problem if critical 
core size and composition are known, then the space-independent solu- 
tion using the buckling may be used. It was found, however, that 
this was not true for subcritical assemblies, i.e. , the space and 
energy dependence of the flux are not separable. 

Further work with larger numbers of energy groups will have to be 
carried out to establish the program's ability to handle successfully many 
group calculations. The unsuccessful attempt at using sixteen groups 
was due to the failure to attain suitable convergence in a reasonable 
amount of time. A better technique than the point successive over- 
relaxation method is obviously needed when dealing with this number 
of energy groups in order to improve the convergence time. 
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VI. RECOMMENDATIONS FOR FUTURE WORK 



In its present form the program has proven to be successful for 
the purpose for which it was intended. Future researchers might choose 
to expand or modify this program in one or more of the following areas. 

A. GEOMETRY AND SPATIAL DIMENSIONS 

Although the program is written for the cylindrical geometry in two 
spatial dimensions (flux cj> is not a function of azimuth angle © ), 

it need not be so limited. The program could easily be adapted to the 
three dimensional cylindrical geometry as well as the two or three 
dimensional rectangular geometry. The addition of a third dimension in 
space would, of course, require considerably more computer core storage 
and hence the future user might have to limit the physical size and 
number of energy groups in a particular problem. 

Another possible additional feature might be to increase the 
number of regions by one or more to include blanket, reflector and 
shielding regions. 

B. CROSS SECTION DATA 

All the results for this program were obtained using three group 
ANL data. Conceivably future users would want to use the most 
current and widely accepted cross section available to them. For 
this reason the program could be modified to handle multigroup data of 
up to twenty-six groups or higher. Once again the user must bear in 
mind that as the number of energy groups is increased, core storage and 
computation time will also increase and probably not linearly. 

C. INNER ITERATION METHODS 

The point successive overrelaxation iterative method gave 
satisfactory convergence rates for the three group problems used in 
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testing the program. However, more rapid convergence of the inner 
iteration is possible [15] using line rather than point successive 
overrelaxation methods and cyclic Chebyshev semi-iterative methods . 
If the future user desires to use more than three energy groups, he will 
almost certainly have to change the program to fit one of these methods 
of solution. 

D. TIME DEPENDENCE 

On a somewhat more far reaching level, this program might be 
expanded to include the time domain for the purpose of studying the 
dynamics of a fast reactor. This could also include a study of space- 
dependent burnup and reactivity feedback effects. These topics are 
currently under investigation by those working on the development of 
the fast breeder reactor. 
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APPENDIX A 



DERIVATION OF LAPIACIAN OPERATOR IN CYLINDRICAL 
COORDINATES AND FINITE DIFFERENCE APPROXIMATIONS 



In order to numerically solve the set of partial differential 
equations 11-31, a finite difference technique was chosen. In 
cylindrical geometry the Laplacian or V ~ operator is approximated 
by transforming the Cartesian coordinates to polar coordinates. 




Y 



X = r cos © ; y = \r sin © 

r = Cx 2 + y 2 ) ' /2 - * 5 © = t ar, -Y. 



A- 1 
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Taking the partial derivations of r and © with respect to x and y in 



Equation A-l results in 



J_ 2 X 

r x ~ 2 <x 2 +y 2 ; y-z. 



- Co s © 3 




2Y _y 

(x a +y z)y 2 _ r 



Sin 0 , 



©x 



Y/x 2 - _ _ 

1 + (>/x) z r a 



s \ r> 6 
r ■> 



© 



y 



. Vx 

i + c y/x) a 




COS Q 

r 



Now consider a function $ C>~ ) z>©) , wherein r and 0 

functions of x and y through Equation A-l. Taking the first partial 
derivative of 0 with respect to x, y, and z utilizing Equations A 
through A-5 gives 

fiy. ~ fie x Y fie ©>< = ~ $© — -p — ? 

fiy - fir Py+ fi 9 &y ~ -fir S'\T\Q + 0e ~ , 

fil. - fi z Z* — 



A- 2 



A-3 



A-4 



A-5 



are 

-2 



A-6 



A- 7 



A-8 
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The second partial derivations are now obtained as follows through 
Equations A-6 through A-8: 



* 



XX 



= C 3 r CoS© - 3 © 

= 0rr COS 2 © + 

S\Y\ Q COS o 
lr 



Sine 



) C 0 



5 \n e 



-V* 



9*< 



20 



e 



r COS© - 0e 

S'm 2 6 
'e e 

Sin© cose 
r ^ 



5 in e 

r 

2. 0r 6 



) 



A- 9 



0yy — ( 3r 5 i r\ © + 



cose 



3 © 



r 

cos 2 © 



) ( (f)v S i n 0 + 0 e 



cos © 



— 0i-rS'm 2 © + 0r p 

+ 2 0TQ COS - V — -2 0 



-V 






ee 



COS 2 © 



r 



A- 10 



© 



Cos© Sin© 



0 



“ az. 



C 0O - z 



A- 11 



Adding Equations A- 9 , A- 10 and A-ll, the Laplacian in cylindrical 
coordinates becomes 

V 2 0 = 0y-r"'-'4 r -0r + 0ee + 0z.z A- 12 

The partial derivatives are replaced by their respective central difference 

2 2 • 

operators 16 with errors of order h where h is the square of the 
distance between adjacent mesh points (See Figure A-l.) resulting 
in 

V 2 0=-^a^ 0(1T+ l,2 ) e)-Z^ 0,21,0) 

+ v T l 0 c r+i ; z,©) - 0 (r-i) 2 . ,©)] + \ 0 Cr,z,©-0 

A- 13 

-2 0Cr,z,©) + 0(r,-z,©+o3 + { 0 Cr > “20(r,z,©) 

+ 0 Cr,z-\ , e) ^ . 
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Figure A- 1 
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where the symbols used in Equation A- 13 are shown in Figure A-l. 

Salvadori and Baron [ 1 6] , as well as any other good text on 
numerical methods, derives finite difference operators with higher 
order errors , i . e . , more accurate approximations . 

It was found that this operator cannot be used for a mesh point 
lying on the longitudinal axis of the cylinder. At these points 
Equation A- 12 has singularities in the second and third terms. In 
this study circular symmetry was assumed, thus Equation A- 12 reduced 
to 

V 2 0 - 0rr + + ^2- . A- 14 

There still existed a singularity in the second term; however, the 
problem under investigation also possessed symmetry with respect to 
the origin, i.e., 0r = 0 at r = 0 . Smith £9] shows that by 
Maclaurin's expansion 

0r C r) = 0r (o) + r r (o) + r 2 0 r rr (0)+« • • , A _ 15 

however, 0r(O) = 0. Therefore, as r approaches zero in the limit 
the term ~~ 0 r in Equation A- 14 is the second partial derivative 

of 0 at r = 0 plus higher order terms. Equation A- 14 can be written 
for mesh points on the longitudinal axis 

S7^0 = 2.0rr + 0 zz. A- 16 

or in finite difference operator form 

V 2 0 = £ 0(r+i ) 7O - 20(1 0Cr-l > 

A- 17 

+ -^> \ 0Cr>z + O - 2 0CV-.Z) a- 0Cr ,z.— \')~] . 
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APPENDIX B 



DEFINITION OF COEFFICIENTS 



In order to simplify the writing of Equations 11-35, they are first 
expanded in order to utilize the sixteen group data available in reference 

7 . 



A typical equation, for example, group three, is written as follows 



[ Yl Co-J ni ,_ 3 +y 3 ^), J o-/,,y] (uz) 



o = i 






^3^ 



+ [ 52 N J (4,2- 3 + %3 Vcrfi'o'] 0 2 (r,Z) 

v — > 



A 3 ( 2 .") 



3 jQ H ' 3 °"fcr ? 3 
> 



i* { $3 (r+i, z) -2<3 3 (r ) z;-hy? 3 (r-. ) x^ 



A 3 ( 3) 



+ 



2.V~h> 



9^3 Cr+i > -z')~9 ! 3C'^-' 1 z ')3' v T'^ f 03 (r >^ +1 '> -2^ 3 0r,z') 



r c 



+ CT 



0 

n 



+ CT, 



in ) 3 



+ cr 



) 3 



-v cr 



vn , 3 



+ cr 



»« •> 3 



+ CT 



ur> 1 3 



8 



A 3 (4} 

+ ^er, 3 + °"5 3 3 9 ^ 3 ^ r >' z -^ 

) 
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+ E 2 T N ^ ( %3 "0 3 cr ; j, 3 )"} 03 ^z) ■+■ ■ • ' 
_ _ ; 



A 3 < 5 ) 



+ T S ^ (X^i a/,\s )] 0 t 6 (r,z) = s 3 (r >^0 

^ v > 

A JOS') 

This equation is then written in the form 

[A 3 C4) + A 3 C3) ( p)]03 (r,z) - v [A3(3)('^ + I7K 



B-l 



03 Cr+ijZ) + £ A3 (3) C 03 (Y 5 2.- 1) + • • • = A3 (0 0 , 0 ) 2 ) 

B-2 

+ A 3 (z) 0a (r,z) -f A 3 (5) 03 (r,z ) + * • * 



+ A3O8) 0*0,2) + S 3 C r > Z. ) 

The general form of this equation is 

[As (g+0 + Ag(g) ( -^)]09Cr } z) + [ A g (< 3 ) cfc + ~ )] 
0*9 Cv+\,-z') + [ As ( 9 ) ( "v 0 ~ aTh 0 sCr->,‘z') + [A 3 (g) 

( 03 (r ,z+i ) + [ Ag <r g ) ( 0 g C v,z-i) b-3 

- Ag(0 0,O,‘ 2 ) - 4 - ... -v Ag ( 9 - 1 ) 0g-' Cr ) 'Z') 

*Ag(g + 2 ) 0 g C r,i) 4- ... + A g ( 18 ) 0, 6 ( r,Z) 

+ S 3 ( r , z ) . 
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This general equation is condensed further to the form 

A 1 jzfgCr-ij 2) + A 2 00Cr ) ^) + A 3 5259 (r+i ,2) +A 4 0g(r, z-i) 

+ A 5 09 Cr, Z+O = AS) 0 ! (TjZ') + * • • dASg-i 09-1 Cr >~Z.) B _ 4 

+ AF9 03 Cr,zH - • • + AF 18 0|6 C'Tj'z.) *+ SgCr^z.") . 

In energy group one. Equation B-3 does not apply. The reason for this 
is that the terms on the right side of Equation B-3 are not fully applica- 
ble to group one. These terms have the following meaning: the sources 
of all neutrons that enter group g due to down scattering from all higher 
energy groups plus neutrons born into group g due to fissions in all 
higher energy groups are represented by the first set of terms on the 
right hand side of Equation B-3 £ Ag( l) 0 1 +- • • • 

+ A3 ( g-l) 0 9 -! Cr ,20 1 . 

In group one these terms must be zero since there are no higher energy 
groups . 

These apparent inconsistencies were accounted for in the compute 
program, i.e. , these terms were set to zero in the first group equation. 
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APPENDIX C 



FINITE DIFFERENCE MESH POINT DESIGNATION SYSTEM 



In order to take advantage of symmetry in the study and hence 
reduce the computer core storage space needed, the following mesh 
point scheme was adopted: a) all points 0 (l,x) = 0 (3,x) and 

0 (x,l) = 0 (x,3) due to symmetry, b) all extrapolated boundary 

points 0 (L,x) = 0 (M,x) = 0 and c) no mesh points have zero 

index due to computer language incompatability . See Figure C-l for 
further details . 
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Figure C-l 
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APPENDIX D 



GAUSS-SEIDEL ITERATIVE METHOD 



In this technique, which is also called the unextrapolated 
Liebmann method, the point single step iterative method, or the method 
of successive displacements, the latest iterative values are used as 
soon as they become available. 

Consider the solution of Poisson's equation 

V 2 0 = S (x>y ) > D-l 



in a rectangular region as in Figure D-l. Assume the solution is proceed- 
ing column by column from left to right. 
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Figure D-l 



If the point at which the equation is presently being solved is 
0 (i,j), the iteration formula is 

0 ( «>p ~ 4 1 0 + 0 + 

D-2 

+ 0^ Ci>j-i') -v i > j+O -v- h 2- $ ( i, j) ^ 

where the standard five-point central difference approximation is 

„ 2 . 

substituted for the V operator, 

POINT SUCCESSIVE OVERRELAXATION 

This method, also called the extrapolated Leibmann method, can 
be thought of as the sum of the Gauss-Seidel iterate Equation D-2 and 

the nth iterate. It is written as 

, Yl+I , v 1 C t\ + i -v-, 

0 C »-i , -v $ 

+ 0 ^ i j j -O 4" 0^ ( i , \+ l ) + h 2 i , j) ~T| D-3 

where w is the overrelaxation or acceleration factor. 

The optimum value of this factor in general is l < w op t < 2. and 
can be computed from 

z 

W = ;■ ■ ■ ■ D-4 

opt 1+ aJ I - X, 2 

where A ( is the largest eigenvalue of the problem resulting from a 
solution using the method of successive displacements, i.e. , the 
Jacobian method. In practice, however, the eigenvalues are seldom 
known and finding them by solving the problem by successive displace- 
ments can be very difficult due to slow convergence. References 8, 9, 

10, 17 recommend an experimental process for determining w ^ f° r a 
particular problem. In this process various values of w are tried until 
one is found that gives the most rapid convergence. This is then dele: mined 
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to be w Q p t ' Reference 16 suggests two other methods for determining 



w 



opt 



The experimental approach described above was used in this 



study. 

The rapid convergence experienced with this method on many 
different problems has proved far superior to the ordinary Gauss-Seidel 
or successive displacements method without overrelaxation. References 
8, 17 contain some examples. 

References 8, 9, 10, 17 relate more explicit details on the 
method of successive overrelaxation in varying degrees. 
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APPENDIX E 



INPUT DATA SYMBOLS 



SYMBOL 


NAME & DESCRIPTION 


UNITS 


REMARKS 


GS 


MESH SPACING 


CM. 


VALUES FROM .75 
TO 2 . 0 HAVE BEEN 
USED SUCCESS- 
FULLY 


GSSQ 


VALUE OF MESH 
SPACING SQUARED 


CM. 2 




C 


TOTAL NUMBER OF 
ELEMENTS OR 
ISOTOPES IN 
THE ASSEMBLY 




UP TO FIVE HAVE 
BEEN USED 


R 


NUMBER OF MESH 
POINTS IN THE 
RADIAL DIRECTION 
PLUS ONE 




ADDITIONAL POINT 
LIES OUTSIDE 
REGION OF SOLU- 
TION BUT IS 
NEEDED FOR DERIV- 
ATIVE BOUNDARY 
CONDITION 


Z 


NUMBER OF MESH 

POINTS IN THE 

AXIAL DIRECTION 
PLUS ONE 




IS IN EFFECT ONLY 
HALF THE 
AXIAL DIMENSION 


SS(I) 


ENERGY 

DEPENDENT 

SOURCE 


r 

NEU ./CM t 


SPACE DEPENDENCE 
MUST BE ACCOUNT- 
ED FOR IN 
DIFFUSION 
EQUATIONS 


S(I) 


GROUP RELATIVE 
MAGNITUDE FACTOR 






DELE(I) 


ENERGY GROUP 
WIDTH 


Mev 




H(I) 


ATOMIC ABUNDENCE 
OF EACH ELEMENT 
OR ISOTOPE 


ATOM/ 

CM. 3 




SIGIN 
(I , J , K) 


INELASTIC SCATTER- 
ING CROSS SECTION 
OF THE Ith ELEMENT 
FROM GP.J TO GP.K 


barns 
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SYMBOL 


NAME & DESCRIPTION 


UNITS 


REMARKS 


SIGTR(IJ) 


THE TRANSPORT CROSS 
SECTION OF THE Ith 
ELEMENT OR ISOTOPE 
AT THE Jth ENERGY 
LEVEL 


barns 




SIGF(I,J) 


THE FISSION CROSS 
SECTION OF THE Ith 
FISSIONABLE ELEMENT 
OR ISOTOPE AT THE 
Jth ENERGY LEVEL 


barns 




SIGC(I J) 


THE CAPTURE CROSS 
SECTION OF THE Ith 
ELEMENT OR ISOTOPE 
AT THE Jth ENERGY 
LEVEL 


barns 




CHI (I) 


NORMALIZED FISSION 
SOURCE FOR THE Ith 
ENERGY LEVEL 






NU(IJ) 


AVERAGE NUMBER OF 
NEUTRONS RELEASED 
PER FISSION OF THE 
Ith FISSIONABLE 
ELEMENT OR ISOTOPE 
AT THE Jth ENERGY 
LEVEL 
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INPUT DATA 
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INITIAL FLUX ESTIMATE 
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